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The modelling of the Milky Way galaxy is an integral step in the study of 
galactic dynamics; this is because knowledge of model parameters that define 
the Milky Way directly influences our understanding of the evolution of our 
galaxy. Since the nature of the Galaxy's phase space in the neighbourhood 
of the Sun is affected by various Milky Way features, measurements of phase 
space coordinates of individual stars that live in this neighbourhood, will bear 
information about such features. Thus, inversion of such measurements can 
^Sj , help us learn the parameters that describe such Milky Way features. 

In the past this has been attempted via the "calibration method", a data 
and computational cost intensive method that is limited by logistical con- 
siderations. We argue in this paper that the calibration method is unsatisfac- 
tory because of reasons both methodological and computational. Here we 
develop a Bayesian inverse problem approach, where we model the avail- 
^ . able stellar velocity information matrix as an unknown function of the Milky 

C/3 ' Way model parameters, where this function is inverted using Bayesian tech- 

niques to predict the model parameters. This unknown function turns out to 

be matrix-variate, which we model as a matrix-variate Gaussian Process. We 

►^ develop a general method to perform inverse, nonparametric learning, us- 

r^ . ing such Gaussian Processes. For the inference we use the recently advanced 

\^ • Transformation-based Markov Chain Monte Carlo (TMCMC). 

^\ I Application of our method to observed stellar velocity data results in esti- 

lO ' mates that are consistent with those in astrophysical literature. That we could 

—^ ' obtain these results using far smaller data sets compared to those required for 

(""^ , the calibration approach, is encouraging in terms of projected applications to 

fT^ ■ external galaxies. 



1. Introduction. Curiosity about the nature of the phase space that we earthUngs hve in, is only 

natural. On astronomical length scales, this translates to the pursuit of the topology of phase space in 

C^ ', the Solar neighbourhood. By phase space of the stars in the Milky Way, is implied W, the space of all 

states that the stars in our galaxy can attain. This reduces to the spatial X and velocity P vectors of all 

the stai^s in the Galaxy, where the spatial location of the a-th star is Xq = (X|° , X2 , X3 )'^, and its 

velocity vector P^ = (Pi''^ P^''\ P^""^ f . 

The phase space of the MiUcy Way in the vicinity of the Sun, has been studied and astrophysical 
modelUng indicates that in this local patch, the phase space structure is sculpted by different dynamical 
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features of the Milky Way (Antoja et al., 2009; Chakrabarty, 2007; Minchev et al., 2009), such as a 
central, rotating, stellar bar in the Galaxy as well as the spiral pattern of the Galaxy. Thus, the phase 
portrait of the solar neighbourhood is affected by the parameters that define these Galactic features. 
Relevant parameters include the orientation of the observer at Eaith (i.e. approximately at the Sun) with 
respect to an axis of the bar, the distance from us (i.e. the Sun) to the centre of the Galaxy, structural 
parameters of the bar and the spiral pattern, as well as parameters that describe the dynamics of these 
features, such as the frequencies of rotation of the spiral pattern il^ and that of the bar Q,b, etc. In fact, 
there are further features that affect the Galactic phase space structure so that the feature parameter 
space demands expansion to include parameters that have hitherto been ignored from the modelling of 
the Galactic phase space. 

Given that these feature parameters affect the local patch of our Galaxy, if phase space coordinates of 
a sample of stars that live in the local patch are measured, such data will bear information about the hand- 
iwork of these relevant features. It then follows that inversion of such phase space information will allow 
for the learning of the unknown feature parameters. This approach has been adopted in the modelling of 
our galaxy, to result in the estimation of the angular separation of the Sun from a chosen axis of the bar 
and the distance of the Sun from the Galactic centre (Dehnen, 2000; Fux, 2001; Minchev et al., 2010; 
Simone, Wu and Tremaine, 2004). The other relevant feature parameters are typically held constant in 
such modelling. 

The modeUing of the solar neighbourhood by Chakrabarty (2007) is similar in that the two-dimensional 
solar location is estimated. The method used in this work is an example of the "calibration method". The 
objective in this work was to predict the solar location given the available stellar velocity measurements 
and simulated phase space data of a sample of stars generated at n different chosen values of the Milky 
Way feature parameter vector S. Each of the n simulated velocity data was employed to generate n 
density estimates in the space of the velocity vectors, via kernel density estimation techniques. The esti- 
mated density functions appeared highly multi-modal and heterogeneous. A non-parametric frequentist 
test was designed to test for the null that the observed data is distributed according to the i-th density 
function that is estimated from the i-th simulated data set that was generated at the i-th chosen value 
of the Milky Way model parameter vector S, (i = 1, . . . , n), (Chakrabarty, 2011). The p- value of the 
used test statistic was recorded for each i. Those choices of the value of S that corresponded to high p- 
values, were considered better supported by the observed data. Hence the empirical distribution of these 
p-values in the space of S, obtained from the n repeated testing, was used to provide intei"val estimates 
of the Milky Way parameter vector. 

Since the best match is sought over a very large collection of training data points - obtained by con- 
structing a grid over the relevant sample space - the method requires computational effort and resources. 
This shortcoming had compelled Chakr-abarty (2007) to resort to an unsatisfactory coarse gridding of 
the space of S. This problem gets acute enough for the method to be rendered useless when the di- 
mensionality of the vector S that we hope to learn, increases. Moreover, the method of quantification 
of uncertainty of the estimate of the location is also unsatisfactory, dependent crucially on the binning 
details of the space of S, which in turn is bounded by cost and memory considerations. 

Thus, an improved method of learning the unknown Milky Way model parameter vector is sought. 
In this paper, we present such a method that allows for the Bayesian non-parametric learning of an 
unknown model parameter vector, given matrix-variate data. Such data is then expressed as a matrix- 
variate function of the model parameter vector where this unknown function is modelled as a reali- 
sation of a Gaussian Process that we discuss here. We demonstrate the effectiveness of this method 
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even with much smaller data sets than were used in astronomical modelling in the past (Chakrabarty, 
2007). That the new method is able to work with smaller phase space data sets, is an important ben- 
efit, particularly in extending the application to galaxies other than our own, in which small numbers 
of individual stars are going to be tracked in the very near future for their velocities, under obsei^va- 
tional programmes such as GAIA (Kucinskas, Lindegren and Vansevicius, 2005; Lindegren et al., 2007, 
http : //www. rssd. esa . int/index .php?pro ject=GAIA&page=index) and PANStarrs 
(Johnston, Bullock and Strauss, 2009); the sample sizes of measured stellar velocity vectors in these pro- 
grammes will be much smaller in external galaxies than what has been possible in our own. The other 
major advantage of this presented method is that it readily allows for the expansion of dimensionality of 
the feature parameter vector of the studied galaxy and is capable of taking enors in the measurements 
into account, within a Bayesian approach. 

Of course, it is of aesthetic interest to learn our location inside our galaxy, but the relevance of this 
exercise goes much beyond. Any modelling of the Milky Way invokes the solar location - in that sense, 
the solar location is of vital importance in the study of the Milky Way. Thus, learning the distance of the 
Sun from the galactic centre, Rq, will crucially constrain the rotational frequency Q,b of the bar in the 
Milky Way, as we discuss later in Section 5.4. The value of ^b affects the understanding of some crucial 
physical processes in the Milky Way. Also, the direction $ = is chosen in our models to coincide with 
the long axis of the central stellai^ bai^ in the Galaxy. Thus, learning the angular separation of the Sun 
from this chosen axis, <I>0, would allow for an understanding of the separation of the line joining the Sun 
to the Galactic centre, from the long axis of the bar; this in turn will allow for improved understanding 
of the effect of the bar on the dynamics of the solar neighbourhood (Englmaier and Gerhard, 1999; Fux, 
2001; Minchevetal., 2010). 

The rest of the paper is stmctured as follows. In Section 2 we discuss the structure of the available 
simulated and real data. The salient aspects of our application are delineated in Section 3 which also 
includes the motivation for a Gaussian Process-based estimation. In Section 4 we present the details 
of the modelling strategy that we adopt. In this context, the relationship between the infomiation and 
the model pai^ameter is presented as a realisation of a Gaussian Process that we discuss. The achieve- 
ment of the posterior probability density of the unknowns, given all data, is arrived at in Section 4.4, 
following discussions about the likelihood and priors on the unknown parameters in Section 4. 1 and Sec- 
tion 8 respectively. The treatment of measrement eixors within the modelling is discussed in Section 4.5. 
Section 5 deals with the implementation of the method to solve the real-life astronomical problem. In 
particular, details of the inference strategem used in our work are given in Section 5.1, Section 5.2 and 
Section 5.3. and Section 5.4 contains results obtained from using the available data. We compare the 
obtained results with the estimates available in the astronomical literature in Section 5.5 and focus on 
the very useful advantage of our methodology (that it can work even with much smaller data sets than 
other available) in Section 5.6. The benefits of our method are delineated subsequently in Section 5.7. 
The paper is rounded up with a discussions section which includes a discourse on other possible areas 
of application of our method. 

2. Data. Essentially, we are attempting an inverse problem, namely that of estimating the unknown 
Milky Way model parameter vector S from the sampled discrete phase space stellar data. The method is 
designed to work with observed real data as well as simulated data, given that in addition to the observed 
data, n simulated stellar phase space data sets are available (Chakrabarty, 2007) at ?i chosen values of 
S, for 4 distinct astrophysical models of the Galaxy. The method is designed to estimate the Galactic 
parameters S G M.'^, that affect the structure of the phase space around the Sun. In our application, we 



4 CHAKRABARTY, BISWAS AND BHATTACHARYA 

restrict the dimensionality of S to 2, i.e. we estimate the coordinates of the radial location Rq of the 
Sun and the solar azimuthal coordinate <I>q. 

It is to be noted that we approximate the geometry of the Milky Way disc as a 2-dimensional disc. 
Thus, we confine analysis to such a two-dimensional spatial geometry, rendering the spatial vector of 
the a-th star {Xf' , X2 )^, and the velocity vector (P|" , P2 )^- Given the disky spatial geometry, we 
find it sometimes useful to use the polar, rather than the Cartesian coordinate system to scan the Galaxy. 
Then the radial coordinate of a particle is R and azimuthal coordinate is $, such that Xi = Rcos $, 
X2 = Rsin^, i? > 0, <^ G [0,360°]. The radial and transverse components of the velocity vector are 
then U and V. 

While the Milky Way model parameters are sought with respect to the centre of the Galaxy, the 
measurements that are available, have been taken by us, i.e. with respect to the Sun. In other words, 
the measured stellar velocities {U^"'\ V^"'^)'^ , a = 1, ... ,j, are with respect to the solar velocity while 
the measured spatial locations of the stars are with respect to the location of the Sun. Furthermore, the 
spatial locations of the sampled stars are within a small distance e from the Sun (Fux, 2001). We interpret 
these data to have been sampled from a circle with radius e centred at the Sun, living in the space of the 
Milky Way disk. Also, the angular- distribution of the sampled stars, inside this 2-ball, is approximately 
uniform. Thus, the summary of the distribution of the measured spatial locations of the sampled stars 
will always coincide with the origin of the used spatial coordinate system, irrespective of what the 
galactocentric spatial location of this origin is, i.e. irrespective of the observer's galactocentric spatial 
location. Thus, these measured spatial locations cannot constrain the sought galactocentric location of 
the Sun. 

So the measurements (from the Sun) of spatial locations of the sampled stars do not bear any infor- 
mation about the unknown (i?0,<l>Q)^, though the measured heliocentric stellar velocities are indeed 
affected by the choice of {Rq, ^0)^- This is because, a given star, if observed from different spatial 
locations in the Galaxy, would appear to move in different ways. For example, if a star appears to have 
a velocity vector directed along the line that joins itself to the obsei^ver at point A on the Milky Way 
disc, this observer will register its velocity to be entirely radial, with zero angular component of the 
velocity vector. On the contrary, had the observer been at a different point B, the velocity vector of this 
star would have registered to have had a radial as well as an angular component, in general. Thus, the 
observed stellar velocities will bear information about the location of the observer, i.e. the Sun. Then 
the available velocity information V can be considered to bear the signature of Milky Way parameters 
S. In principle, beyond just the galactocentric solar location, if there are Milky Way parameters that 
physically affect the structure of the phase space, velocity information drawn from that phase space will 
bear the signature of such Milky Way parameters. 

In light of the above discussion, the measured data used in our work is considered to comprise stel- 
lar velocity information, {{Ua,ya)'^Va=i^ where a total of j number of stars are observed. Thus, the 
measured velocity information is a matrix but in our modelling, we find it more useful to consider the 
j number of measurements of the /c-dimensional velocity vectors, as a j/c-dimensional velocity vector 
v(*"=**). In our application to learn the Milky Way model parameters, the dimensionality of the measured 
velocity vectors is 2, i.e. k=2. 

As mentioned above, we supplement information about the unknown Milky Way parameters with 
simulated information that is generated from dynamical simulations of 4 astronomical models of our 
Galaxy, at n distinct choices of values of S. In fact, the simulated data generated for the q-ih as- 
tronomical model, corresponding to chosen values, s^, . . . , s* of 5, comprise the training data Vs, 
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g = 1, 2, 3, 4. At the i-th chosen value, s* of S, let there be j number of fc-dimensional velocity vectors 
simulated, i = 1, . . . ,n. Thus, we will expect such simulated information to be in the form of a series 
of n matrices, each of which is j x fc-dimensional. However, as with the real data, we find it useful in 
our modelling to treat each simulated velocity matrix as a jfc-dimensional vector, and there are n such 
generated velocity vectors, generated at the n different choices of the model parameter vector S. From 
this viewpoint, the entire simulated velocity information is a matrix Vs, with n rows and jk columns, 

( n X i k^ 

i.e. the simulated velocity information can be considered to be the matrix Vl , such that the i-th row 
is the jfc-dimensional simulated velocity vector Vj, generated at the i-th chosen value s* of S, where S 
is itself a d-dimensional vector. 

In our work, the simulated Vg is referred to as the training data while the real or observed data v^*'^**) 
is referred to as test data. This is because the simulated data are generated at known values of S while 
the model parameter supported by the real data is unknown. The set {s*, . . . , s*} of the n chosen values 
of S at which the simulated velocities are generated, is refen^ed to as the design set and Sj as the i-th 
design vector. 

The inverse method that we advance below is a generic methodology that can be applied to estimate 
the d-dimensional vector S, where d is bounded from below by the number of parameters in the astro - 
physical models that can be called upon to bear influence on the phase space structure of the Galaxy 
in the solar neighbourhood. As we have said above, we demonstrate an application of our developed 
learning technique to the learning of the 2-dimensional solar location vector. The reason for restricting 
our application to the case of d=2 is the existence of stellar phase space data generated in dynamical 
simulations of astrophysical models of the Milky Way which involve scanning over chosen guesses for 
Rq and $0, with all other feature parameters held constant. If simulated data distinguished by choices 
of other Milky Way model feature parameters become available, then the implementation of such data as 
training data will be possible, allowing then for the learning of Milky way model parameters in addition 
to Rq and <1>q. The methodology that we advance is indeed a generic one, with computational costs 
being the only concern in extending to cases of d > 2; that extending to a higher dimensional S only 
linearly scales computational costs, is discussed later in Section 6. 

Actually, the learning of the unknown S is attempted 4 times, each time using the observed data v^**^^*) 
and one of the 4 distinct Vg' data that are generated at the chosen design set, using the 4 astronomical 
models that were considered by Chakrabarty (2007), g' = 1, 2, 3, 4. Above, the training data Vs' is the 
n X j /c-dimensional matrix that is simulated using the q-th astrophysical model of the Milky Way, such 
that its i-th row is the i-th velocity vector Vj that is generated at the choice of S = s* i = 1, . . . ,n. 
The choice of the base astrophysical model is distinguished by the ratio of the rotational frequencies of 
the spiral to the bar, Qs/^b- That this ratio is relevant to the phase space structure of the model of the 
Galaxy is due to the fact that 0,s/^b can crucially control the degree of chaos in the Galactic model. For 
example, it is well known in chaos theory that when 0,s/^b is such that one of the radii at which the 
bar and the stellar disk resonate, concurs with a radius at which the spiral and the stellar disk resonate, 
global chaos is set up in the system (G. Walker and J. Ford, 1969). Chaki^abarty and Sideris (2008) have 
corroborated that the degree of chaos is maximal in the astrophysical Galactic model marked by such 
a ratio (r2s/il6=22/55). They report that in models marked by slightly lower (Jig /0;,= 18/55) or higher 
(Us/^b = 25/55) values of this ratio, chaos is still substantial. In the Galactic model that precludes 
the spiral however, chaos was quantified to be minimal. It is these 4 states of chaos - driven by the 4 
values of ^s/^b - that mark the 4 astrophysical models as distinct. Thus, the 4 astrophysical models 
lead to model phase space structures that are differently chaotic. This results in 4 distinct simulated 
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velocity information sets Vs , Vs , Vs , Vs that bear the effects of such varying degrees of chaos, 
each generated at the chosen design set {s^, . . . , s*}. As attempts to learn the unknown Galactic model 
parameter vector are undertaken with one astrophysical base model at a time, each such attempt is 
independent of the other. Consequently, we simplify the notation in our discussion of the method below, 
by dropping the reference to the index q on the training data sets hereafter. 

Repeated learning of the unknown s("'^"') using the real measurements v*-**^**^ and one of these sim- 
ulated data sets will allow for the examination of the influence of a chaotic phase space structure on 
such learning. For instance, if the implementation of simulated data sets that are generated from base 
astrophysical models marked by similar chaoticities, result in learnt values of ^("sw) j^at are widely 
disparate, that would be cause for concern. At the same time, a method that is unable to perform suc- 
cessful learning of s^"-^^^ in the face of extreme chaos - such as the method of Chakrabaity (2007) - is 
viewed as weak. We show below that in our method, the robust model uncertainty estimation within a 
Bayesian framework allows for the learning even in the face of global chaos in the base astrophysical 
model. Details of the dynamical simulations performed on the 4 astrophysical models are given in the 
supplementary section S-1. 

3. Case study. In this work, our interest lies in predicting the value of the unknown Milky Way 
model parameter s("'^"') that supports the real data vector v^**^**), given the training data matrix Vs that 
is generated at chosen values s^, . . . , s*, of the Milky Way model parameter vector S. We define the 
chosen design set as {s|, . . . , s*}. The training data is employed to train the model for the unknown 
functional relationship between the available velocity information and the set of values of the model 
parameter vector that support this information. We present the learning of the solar position in the Milky 
Way disk by inverting real stellar velocity measurements in conjunction with simulated stellar velocity 
data generated at a chosen design set. Then we are interested in the inverse problem of predicting s("^"') 
given Vs, v^*^''*) and the trained model for the aforementioned functional relationship. 

For this inversion, one possible way is to first estimate this function from the data, using well-known 
nonparametric techniques, and then use numerical methods to achieve the inverse solution. The com- 
plexity of both these exercises - particularly the latter - increases as the dimensionality of the function 
space increases, unless simplifying assumptions are made. In fact, modelling high-dimensional func- 
tions using splines/wavelets require restrictive assumptions and in particular, may fail to adequately take 
into account the coiTclation structure between the component functions. Furthermore, this approach has 
been criticised on the ground that it ignores parameter uncertainty, in that the approach provides only a 
point estimate of the set of model parameters. 

A Bayesian equivalent of this approach could involve the modelling of the unknown underlying func- 
tion with splines or wavelets, and then computing the posterior distribution of s^"''^^', conditional on the 
training and the test data sets, after marginalising over the other model parameters. This computation, 
as well as the marginalisation, are in principle straightforward to implement using well-known Markov 
Chain Monte Carlo (MCMC) methods. Thus, the Bayesian approach avoids the difficulties of nonpara- 
metric model parameter estimation and subsequent inversion. Additionally, the Bayesian approach does 
not ignore parameter uncertainty while predicting the unknown location. Furthermore, priors on the un- 
known parameters can bring in extra information into the model. Thus, in principle, a training data set of 
even smaller size than that required in the classical approach, will be adequate in the Bayesian paradigm. 
However, as already discussed above, modelling high-dimensional functions using splines/wavelets is 
somewhat restrictive in that this approach fails to provide an appropriate coixelation structure for the 
component-wise functions. 
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In this paper, we circumvent tliis problem by introducing a Gaussian Process approach to modelling 
the high-dimensional functional relationship between available information vectors V and the value of 
S that supports this information. The velocity information that results from collating such information 
vectors is matrix-variate. The functional relationship between a general information matrix V and the 
corresponding matrix of values of S is modelled in our work with a Gaussian Process. In the face of the 
matrix-variate nature of the information, this functional relationship needs to be modelled as a realisation 
from a matrix-variate Gaussian Process (GV)- 

Indeed, had we been treating the infomiation V generated at a value of S^ as a j x /c-dimensional 
matrix, we would have modelled the functional relation between S and V by a j x /c-dimensional 
matrix- variate QV. In that case, one realisation from such a matrix-variate QV would have been the n 
number of j x /c-dimensional matrices containing velocity values that are generated at n values of S. 
However, as said before, in our treatment, the j number of /c-dimensional velocity vectors generated at 
any value of S is treated as a j /c-dimensional vector. In this paradigm, each of the n such infomiation 
vectors that will be considered to be a realisation from a vector-variate QV of dimensionality jk. Thus, 
our view allows for the reduction of the problem from a matrix-variate QV to a vector-variate QV. 

So we treat the general j x /c-dimensional stellar velocity information as a j /c-dimensional information 
vector V'--''^^^). Then we model V as a function of the value s of a general model parameter vector S at 
which such information is generated. Thus, we write V = ^(s) where s is the value of S that supports 
this information vector V. To be precise, in our application, each stellar velocity vector is 2-dimensional, 
i.e. /c=2. Also, we seek to learn a model parameter vector that is 2-dimensional, i.e. (i=2. In the Bayesian 
approach that we adopt, a much smaller j (=50) allows for inference on the unknown value s^"'^"'^ of the 
Milky Way model parameter vector, than j ^^3000 that is demanded by the aforementioned calibration 
approach used by Chakrabaity (2007). The simulated data presented in Chakrabaity (2007) that we use 
here, is generated at 216 distinct values of S, i.e. n=216. 

Thus, our design set comprises 216 chosen values s^, . . . , s* of the model parameter vector. For each 
of the 4 base astrophysical models, at each chosen s*, 50 2-dimensional stellar velocity vectors are 
randomly drawn from the velocity data that are generated from dynamical simulations performed at that 
value of S. Here, i = 1, . . . , 216. As explained before, the 216 pieces of velocity information generated 
in this way are treated in our work as vectors vi , . . . , V216, where Vj has dimensionality 50 x 2 = 100 for 
each i. Then the training data in our work comprises all such vectors and is represented as Ps ^ ' . 
The real data is treated in our work as the 100-dimensional vector v^*^**). The i-th design vector is s^ 
which is a (i-dimensional vector V i = 1, . . . , n. 

It will be seen in Section 5.4 that, in spite of our relatively small number of stars in the data sets, 
the summary of the posterior distributions of the unknown solar position is consistent with existing 
knowledge in the astrophysics literature. The other very important contribution of our work is that it 
allows for the learning of a high-dimensional model parameter vector. 

4. Modelling. In line with the discussion in Neal (1998), in our application we have n j A;-dimensional 
simulated velocity vectors (inputs), each generated at a chosen value of the model parameter vector (tar- 
get), i.e. we have the n observations (vi, s\), . . . , (v„, s*), (where s* is a known d-dimensional vector, 
\f i = 1, . . . ,n). Our aim is to predict the unknown Milky Way model parameter vector s^"^"') when 
we have only information about the input in the form of the real velocity vector v^*^'^*). A predictive 
distribution of gC"*^™) is sought, conditioned on the real velocity vector and the training data generated at 
the known design vector values, namely at s*, . . . , s*. This is achieved in an approach in which we treat 
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the input - that is to begin with matrix-variate - alternatively as a vector of corresponding dimensions, 
allowing us to then consider the functional relation between the inputs and the targets as a realisation 
from a vector-variate QV. In this section we ignore measurement errors and present our model of these 
n vector-variate functions. Later in Section 4.5, we delineate the method used to take measurement 
uncertainties on board. 

We define the full information matrix V^"^-''^^ as one with the i-th row given by the velocity infor- 
mation vector (Vj ^ Y' , where i = 1, . . . , n. Then Vj is supported by the value Sj of the model 
parameter vector S, so that the equation V = ^(s) expresses the functional relationship between the 
two variables V and S. We express the i^-th element of the infomiation matrix V as Vu which is equal to 
the i-th component function ^^(sj) of the jfc-dimensional $,{si). Then S,e{si) is a scalar valued function 
of the vector Sj, where i = 1, . . . , n and i = 1, . . . , jk. In other words, 



y^x V22 • • • V2 jk 

\ Vnl K2 • • • VnjkJ 



6(^2) ^2(52) ••• ijk{s2) 
\ il{Sn) i2{Sn) ••• ijk{Sn) / 



Thus, the i-th row of the V^"^-''^^ matrix is the jfc-dimensional vector Vj that is set equal to ^(sj) := 
(^i(sj), . . . , ^jfc(sj))^' V i = 1, . . . , 71. We model the jfc-dimensional function ^(•) as a jTc-dimensional 
vector-variate QV. One realisation from this QV is {^(si), ^(52), . . . , ^(s„)}. Then the joint distribu- 
tion of {^(si),^(s2), . . . ,^(sn)} is matrix normal, with adequate parametrisation. We represent this 
as 

(4.1) {^(^l), ^(^2), . . . , ^{Sn)} ~ MMnjkit^, A, fi), 

where the mean matrix of this matrix nonnal distribution is the n x jfc-dimensional matrix /i, the left 
covariance matrix is the n x n-dimensional A and the right covariance matrix is the jkx jA:-dimensional 
matrix fi. Before describing these individual matrix-variate parameters of this distribution below, we 
note that Equation 4.1 is the same as saying that the likeUhood is matrix normal. 

In our work we choose to implement the popularly used square exponential covariance function 
(Rasmussen and Williams, 2006; Santner, Williams and Notz, 2003; Scholkopf and Smola, 2002). This 
covariance function is easy to implement. It renders the sampled functions smooth and infinitely dif- 
ferentiable while the covariance structure itself is dependent only on the Euclidean distance between 
two values of the function variable so that the covariance is invariant to rotations and translations in 
the space of this variable. It is hard to preempt if such a stationary and isotropic covariance function 
is too restrictive to properly model correlations for the application at hand. While acknowledging such 
a possibility, we resort to relaxing the choice of a zero mean function though that is a popular choice. 
Instead we choose to define the mean function in a way that is equivalent to the suggestion that the 
velocity information is viewed as centred around a linear model with the residuals chai^acterised by a 
vector-variate QV (A. O'Hagan, 1978; Cressie, 1993). We then integrate over all such possible global 
intercepts to arrive at a result that is more general than if the mean is fixed at zero. An advantage of the 
non-zero mean function is that in the limit of the smoothness parameters (characterising the smoothness 
of the functions sampled from this QV) approaching large values, the random function reduces to a lin- 
ear regression model, which is plausible. In models with a zero mean function though, in this limit of 
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very large smoothness, the random function will concur with the errors, which appears less intuitively 
pleasing. 

The non-zero mean function of the QV is defined to be /x(-) := HB, where 

jjT .- [^(™xl)(si),...,^(™xl)(s„)], with 

(4.2) ^(™>^i)(.) = (/ii(.),...,/i™(-)f 

For e = 1, . . . ,m, he{-) could be any function. However, in accordance to the suggestion that such a 
non-zero mean function be expressed in terms of a few basis functions (Rasmussen and Williams, 2006), 
we choose to fix this functional form as h{-) := (1, ■)'^ for all values of S. (A similar construct was used 
by Blight and Off (1975) who performed a ^P-based polynomial regression analysis). We also choose 
to set TTi = d + 1. Thus, in our treatment, /i(-) is a (d + 1) -dimensional matrix. The extent to which this 
chosen mean function deviates from being zero is borne by the coefficient matrix B which is the other 
factor in the definition of the mean matrix. We define B as 

(4.3) S = (/3n,...,/3ji,...,/3ifc,...,/3,fc) 

where for p = 1, . . . ,j, p' = 1, . . . , k, (3 , is an ?n-dimensional column vector. As we choose to set 
m = d + 1, B is a. matrix with d + 1 rows and jk columns. 

The amplitude of the covariance function is characterised by the matrix Q which is jA;x jTc-dimensional 
and is defined as 

(4.4) n = E ® C 

where S is the kx k covaiiance matrix for a given star in a row of the matrix- vari ate function ^(•). Thus 
S corresponds to the covariance amongst components of the velocity vector of each star for a given 
value of S. On the other hand, C is the j x j covariance matrix for a given component of the velocity 
vector for the stars at a fixed value of S. 

The left covariance matrix of the matrix-nonnal distribution defined in Equation 4.1 is A which is 
n X re-dimensional. Given the choice of squai^e exponential covariance function, it is 

^(nxn) ._ [a{-,-)], where 

(4.5) a{s,s') = exp{-{s-s'fQ{s-s')}, 

for any 2 values s and s' of S. Here, Q^'^^'^' represents the inverse of the scale length that underlies 
correlation between functions at any two values of the function variable. In other words, Q is the ma- 
trix of the smoothness parameters. Thus, Q is a matrix that bears information about the smoothness 
of the sampled functions; it is a diagonal matrix consisting of d non-negative smoothness parameters 
denoted by 6i, . . . , 6^^. In other words, we assume the same smoothness for each component function 
of ^(-j. This smoothness is determined by the parameters 6i, . . . , 6^^. We will learn these smoothness 
parameters in our work from the data. Of course, though we say that the smoothness is learnt in the 
data, the underlying effect of the choice of the square exponential covariance function on the smooth- 
ness of the sampled functions is acknowledged. Indeed, as Snelson (2007) states, one concern about 
the square exponential function is that it renders the functions sampled from it as artificially smooth. 
An alternative covariance function, such as the Matern class of covaiiances (Matem, 1986; Snelson, 
2007; Tilmann Gneiting and William Kleiber and Martin Schlather, 2010), could give rise to sampled 
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functions that are much rougher than those obtained using the square exponential covariance function, 
for the same values of the hyperparameters of amplitude and scale that characterise these covariance 
functions(see Chapter 1 of Snelson's thesis). 

Let ujj-ii denote the (r, £)-th element of il, Cri the (r, ^)-th element of C and let art denote the (r, ^)-th 
element of 5] . Let the ^-th component function of ^ ( • ) be ^£ ( • ) with I = mik + m2, where i = I, . . . ,jk 
and 1712 = 1,2, . . . ,k, rrii = 0, 1, . . . , j — 1. Then the correlation between the components of ^(•) yields 
the following correlation structures: 

Cmik+m2{Si),Cm[k+m2i^i)) = ' V ^2,1 and ITli ^ m^ 

(4.7) corr (Cmife+m2(sj)>Cmifc+m^(si)) = ^===^= Vmi,i and m2 / m'2 

(4.8) corr (Cmifc+m2(si)'?mifc+m.2(*i)) = Vz,mi ^ mi,m2 7^ r^a 

(4.9) corr{^(,{si),^g{s2)) = a(si,S2)V £ and si / S2 

4. 1 . Likelihood. Along these lines, the f-th row of the training data Vs that is itself n x j fc-dimensional, 
can be expressed as a function of the i-th chosen design vector s* of the Milky Way model parame- 
ter vector S. Here the i-th row of Vg is the jA;-dimensional vector ^rf that is simulated at S = s*, 
i = I, . . . ,n. Thus, we can write Vj = ^(s*). Our ultimate aim is to learn the unknown value of S at 
which the real data v^*^''*) is generated. In this attempt, we write the probability distribution of the train- 
ing data Vg, given the parameters that define the distribution of ^(•). In order to achieve this likelihood, 
we define 

• the n X jTc-dimensional mean function Hj^B, where the linear form of the mean structure is 
contained in H^^""""^ ■- [/i(™xi)(s*), . . . , /i(™^^)(s;;)] and the coefficient matrix B is defined 
in Equation 4.3. 

• thesquareexponentialfactor in the covariance matrix A)^^" := [exp{ — (s* — s'*)^Q(s* — s'*)}] 
(see Equation 4.5). 

• the amplitude of the covariance matrix Cl^i^^^''^ = JiC^x'^) (^ (^Oxj) (see Equation 4.4), 

Then it follows from the matrix normal distribution of in Equation 4.1, with mean function defined in 
Equation 4.2 and Equation 4.3, and covariance matrix defined using Equation 4.5 and Equation 4.4, 
that Vs is distributed as matrix normal with mean matrix HdB, left covariance matrix Ad and right 
covariance matrix fi, i.e. 

(4.10) [Vs I B, C, S, Q] ~ MMn,jk{HDB, Ad, n) 

Thus, using known ideas about the matrix normal distribution - see Dawid (1981), Carvalho and West 

(2007) - we write 

(4.11) 

[Vs I B, C, S, Q] = -r^-—, exp |-itr [n~\Vs - HdB)^A~^\Vs - HdB)] ] 



The interpretation of the above is that the r-th row of [PslB, S, C, Q] is multivariate normal with 
mean corresponding to row of the mean matrix HdB and with covariance matrix fl. Rows r and i 
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of [PslB, S, C, Q] has covariance matrix a{sr, S£)Q. Similarly, the ^-th column of it is distributed as 
multivariate normal with mean being the i?-th column of HdB and with covariance matrix uji^iAi), 
where Ur/ denotes the (r, £)-th element of fi. The covariance between columns r and £ is given by the 
matrix oj^/^d- 

4.2. Predicting s("<^"'). In order to predict the unknown MiUcy Way model parameter vector s("'^"') 
when the input is the measured real velocity vector v^*^***), we write the posterior predictive distribution 
of s^^'^^\ given the real velocity matrix v^**^**), the training data Vg and all the known targets, i.e. the 
chosen design vectors s*, . . . , s*. This posterior predictive is usually computed by integrating over all 
unknown model parameters. While it is possible to analytically integrate over B and C, I] and Q cannot 
be analytically integrated out. In fact, we find it useful to learn the d smoothing parameters (the diagonal 
elements of Q), given the data. Thus, one useful advantage of our method is that the smoothness of the 
process does not need to be imposed by hand, but can be learnt from the data, if desired. 

Given our learning of s("'^"'), I] and Q, we rephr^ase our motivation as seeking to compute the joint 
posterior probability of s(""="'), Q and S, conditional on the real data and the training data, for a choice 
of the design matrix. In fact, we achieve a closed form expression of this joint posterior of s^"^"'^ Q and 
S, by integrating over the other hyperparameters, namely, the amplitude of the mean function (B) and 
the matrix C that bears information about covariance among stars, for a given component of the velocity, 
at a fixed value of S. From this closed form expression, the marginal posterior probability densities of 
Q, S and any of the d components of the sC"^*") vector can be obtained, using the transformation based 
MCMC sampling method (Dutta and Bhattacharya, 201 1) that we adopt. 

For a given choice s*, . . . , s* of the design vectors, the posterior distribution [s("'^"') , S, Q| v^*"^**) , Vg] 
is sought, by marginalising the posterior [s(**^**), S, Q, B, Clv^*'^**),!?^] over the process matrices B 
andC. 

4.3. Priors used. We use unifomi prior on B and a simple non-informative prior on C, namely, 
7r(C) oc| C |~u+i)/2. As for the priors on the other parameters, we assume uniform prior on Q and use 
the non-informative prior vr(I]) oc| I] |-('«+i)/2 The prior information available in the astrophysical 
literature (Chakr^abarty, 2007) on s("<^"') suggests uniform piiors on all components of the s("'^"') vector 
(see Section 5.4 for greater details). 

4.4. Posterior of s^"''^^' given training and test data. Since our interest lies in the prediction of 
s ("£«") ^ given the real test data and the simulated training data, as well as in learning the smooth- 
ness parameter matrix Q and the matrix Xl that bears the covariance for a given star, we compute 
the posterior distribution [s^"''^^' ,Q,'S \ \'^^'^^^\'Ds]. As expressed above, we achieve this by writing 
[s("e"'), B, C,Q,T.\ v(*^^*), P,] and marginahse over B and C. 

To construct an expression for this posterior distribution, we first collate the training and test data 

to construct the augmented data set Pj„„ = (v^: . . . :v^:(v''*'^'^*))^). Then the set of values of the 
model parameter vector S that supports Vaug is {s*, . . . ,s*, s^"^"')} of which only si"*^"") is un- 
known. We define /i"i5"+^^'''"^ := [/i('"^i)(s|), . . . , h^"'''^\s*J, /i('»xi)(s("e"'))], where our choice 
of the functional form of h{-) has been given in Section 4 and we also set ?n, = d + 1. Also, we 
define Aii" '^^"' '' := [expj — (sj — s'-,)'^Q(s'- — s'-,)\] where s' and s'-, are members of the set 

K,...,<,s("--)}. 

The priors used on B, C, Q, XI and sC"^'") are listed in Section 8. Using these, we observe that 
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[s("'="'),Q,B,S,C| v(*^^*),P,] OC [Vaug \ B, S, C, Q, s^""*")] [B, S, C, Q, s^'^'^"')] 



jfc, j(n + l)+fc + l fc(n + l)+j + l 

A-n 25] 2 C 2 

I J-^auq I II II 



X exp 



X exp 



■;^(r{n-'l)J„,M„„5D„„,} 



4'-{""{ 



\ l-^aug L>aug ^aug J 



Hjy 



'aua^Vaug^'^^a] 



H 



T^aug T>au 



H 



v„ 



B - {Hl^ 



A^ H-r, 



H-r, A^ 2?o,,g > > 



'»a«9^0a„„'^««9 



(4.12) 



where we have recalled Equation 4. 1 1 . As in the definitions of Hjy^^^ and Ax>„„g , the matrix Maug 
is defined at the given choice of the design matrix, augmented by the unknown Milky Way model 
parameter. 

Marginahsing [s("'^"') , Q, S, X), C | v'^*'^'^*) , Vs] (as expressed in Equation 4.12) over B and C yields 
the joint posterior [gC"^'") , Q, S | v^*^***) , Vg] , the form of which is obtained as follows: 



I f[s^''^'"\Q,B,'S,C\ v(*^^*),P,]dSdC 



ik m -I 1k i(n4-l — m)-\-k4-l ^ 

« \Av.J~-\{Hr,^^y{Ar,^^^}-HHT,^J-^ x |S| 2 \^n + I - m)kCGLS., 

(4.13) 



augl 



2 



We thus obtain a closed-form expression of the joint posterior of (s^"'^"'), Q, S), given training and test 
data, for a given choice of the design matrix (Equation 4.13), up to a normalising constant. The QV 
prior is strengthened by the n number of samples taken from it at the training stage. We sample from 
the achieved posterior using MCMC techniques to achieve the mai^ginal posterior probabilities of Q, 
5] or any component of s^"'^^\ given all data. In fact, we implement the transformation-based MCMC 
method recently by Dutta and Bhattacharya (2011). 

4.5. Errors in measurement. In our application, the errors in the measurement of the stellar veloci- 
ties are small and will be ignored for the rest of the analysis. In general, when errors in the measurements 
that comprise the training data and the test data are not negligible, we assume Gaussian stellai^ velocity 
errors ej, in vj, with t = 1, 2, . . ., such that ej ~ Mjk{0, <i), where (^ = Si ® S2; Si, S2 being positive 
definite matrices. If both Si and S2 are chosen to be diagonal matrices, then (^ is a diagonal matrix; as- 
suming same diagonal elements would simplify <;^ to be of the form 93 x /, where I is the jk x jk-th order 
identity matrix. This enor variance matrix <;^ must be added to Q, before proceeding to the subsequent 
calculations. TMCMC can be then be used to update <j. 

Our fully Bayesian approach quantifies the prior information in terms of the Gaussian process and 
coherently processes the information contained in the prior and the data to yield results that are in agree- 
ment with those obtained by previous approaches using very large datasets. Moreover, unlike previous 
approaches, our approach allows for neat probabilistic conclusions a posteriori. 
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5. Analysis of real data. In our application, 5 is a 2-dimensional vector, i.e. d=2. Also, P is a 
two-dimensional vector, i.e. k=-2. We consider 50 simulated/real stellar vectors at any value of S i.e. 
j=50. Lastly, we work with n design vectors where n = 216. 

5.1. Posterior inference using TMCMC. Motivated by the fact that the performance of traditional 
MCMC methods - including the Metropolis-Hastings algorithm - can be less than satisfactory in high 
dimensions, both in teiTns of convergence and computational time, Dutta and Bhattacharya (2011) pro- 
posed a novel methodology that works by constructing proposals that are deterministic bijective trans- 
formations of some arbitrary random vector drawn from a chosen distribution. Dutta and Bhattacharya 
(2011) show that for additive transformations, the TMCMC-based acceptance rate decreases at a slower 
rate compared to block random walk Metropolis algorithms. Furthermore, TMCMC includes the hybrid 
Monte Carlo (HMC) algorithm as a special case and in one-dimensional situations while it boils down 
to the Metropolis-Hastings algorithm with a specialised proposal mechanism. 

In TMCMC, the random vector of which a proposal density is a transformation of, can be chosen 
to be of dimensionality between 1 and the dimensionality of the pai^ameters under the target posterior. 
High-dimensional parameter spaces are explored by constructing bijective deterministic transformations 
of a low-dimensional random vector. Interestingly, whatever the transformation, the acceptance ratio in 
TMCMC does not depend upon the distribution of the chosen random vector. For our purpose, we shall 
consider TMCMC based on additive transformations, since Dutta and Bhattacharya (2011) show that 
these transformations require far less number of "move types" compared to non-additive transforma- 
tions. 

5.2. TMCMC algorithm. TMCMC allows updating the entire block (s^"^"'), Q, S) at the same 
time. The algorithm is as follows. 

(i) Initialise the unknown quantities by fixing arbitrarily initial values (s^"'^^''^\Q^^','S^^'j. In 

our case, s^"'^''"'^^ = (s^^' , . . . , sj^'^^' ), Q^^' is characterised by the initial values of the 

d smoothness parameters, which we denote hy b := (b^ , . . . ,b^ )'^ and "S^^' denotes the initial 
choice of the k x k matrix S. SI is decomposed into LL^, where L is the appropriate lower- 
triangular matrix 

(ii) Let (f = ((s^"'^'"))^, b^ , (L*)'^)'^, where L* denotes the column vector consisting of the non-zero 
elements of L. 

(iii) Next we propose e ~ 5'(")-^{e>o}' where g{-) is some arbitrary distribution, and / denotes the 
indicator function. In our applications, we shall choose g{-) = A^(0, 1), so that, e > is drawn 
from a truncated noniial distribution. 

(iv) Assume that at iteration t, the state of the unknown parameters is (gl"'^"''*)^ Q^*\ S'-*^) := if^^\ 
Update yj^*) by setting, with probabilities ttj and (1 — ttj), ip- = tp- ± CjC (forward transfor- 
mation) and (f- = (f- —CjE (backward transformation), respectively, where, for j = 1, . . . ,d, 
TTj are appropriately chosen probabilities and Cj are appropriately chosen scaling factors. Assume 
that for ji G U, if, gets the positive transformation, while for J2 G W^, ip, gets the backward 
transformation. Here U UW = {I, . . . ,d*}, where d* = 2d + !<!s+ll_ 
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(v) We accept the new proposal yj(*+^) with acceptance probability 



■<fi — i^ii^i ^ i) rf Ff r; 7 ^ < v 



(5.1) q;<o = mm < 1, -^ fF 7T 7 ^ ^' 



"]2) ) 

aug)\ '^ 



X 



where r^ denotes the ratio of [| AD„„^(sa„g)| '2 |/i-|;^^^(sang)A^^^^(sa„g)fl"D„„s(s 
|5]| 2 |(n + 1 — ^Ji)kCGLS,aug\ ^ , evaluated at the new value {i^^^^^') and 

the current value (sp^^^) of (^ respectively. We only need to bear in mind that the acceptance prob- 
ability is zero if 6j < for any j or if any diagonal element of L is negative. 

5.3. Details of our implementation ofTMCMC. In our work we use the following proposal mecha- 
nism: at each iteration we generated e ^ A/'(0, l)/{e>o}. Then, with equal probabilities, for ^ = 1,2, we 

set s^^^' = s^^^' ± e and 6^ = 6^ ± O.OOTe. For updating 5], we decompose S = LL' , but 
here we set the scale factor to be 0.07 for all the non-zero elements, using the additive transformation 
with equal move-type probabilities. 

Following Section 5.2 we update all the unknowns in a single block, using additive transformations of 
a one-dimensional random variable. This, along with adequate choice of the scale parameters, ensures 
fast computations and reasonable convergence. 

For proper choices of the scale parameters of the additive transformation and the initial values of 
the parameters we conducted several initial "pilot" TMCMC runs of length around 100,000, starting 
with arbitrary initial values and guesses of the scale parameters such that all the runs converged to the 
same distribution as indicated by informal diagnostics such as trace plots. For the final TMCMC run, 
we chose those scale parameters that yielded the best convergence (with respect to empirical diagnostics 
such as trace plots) among the pilot runs, and selected the final values of the parameters obtained in 
this best pilot run as the initial values for the final run of TMCMC. The pilot runs yielded the proposal 
mechanism that we worked with. 

Once the proposal mechanism and the initial values are thus decided, we then discarded the next 
100,000 iterations of our final TMCMC run as burn-in and stored the next 1,000,000 iterations for 
inference. For each model it took approximately 6 hours on a laptop to generate 1,100,000 TMCMC 
iterations. 

5.4. Results from real data. The prior information available in the astrophysical literature (Chakrabarty , 
2007) for the radial distance of the Sun from the galactic centre, Rq G [1.7, 2.3] and the angular sepa- 
ration of the Sun-Galactic centre line from the chosen long-axis of the elongated bar, $0 e [0°, 90°]. 

The radial coordinate values are expressed in units of the "co-rotation radius" Rcr of the central bar 
in the Milky Way disc. This is the radius at which the rotational frequency Q.i^ of the rigidly rotating bar, 
equals the radius-dependent rotational frequency Q.{R) of the stars at distance R from the centre of the 
Galaxy, which in turn is determined by the choice of the model for the gravitational field of the Milky 
Way disc that is one of the inputs in all 4 astrophysical models, the dynamical simulations of which 
produce the 4 training data sets that we implement in our work. In the model for the gravitational field 

of the Milky Way disk that we work with, this rotational frequency i7(i?) of the stars is defined as — , 

R 
where vq is a constant that is set to unity in the simulation models of Chakrabarty (2007) and ilf, is set 

to unity; thus Rcr = 1 in this scaled simulation model. All spatial lengths in our work are expressed in 

units of Rcr- 
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Astrophysical literature (Binney and Tremaine, 1987) offers the prior infonnation that the radial lo- 
cation of the Sun is around URcr, which corresponds to the location of another major resonance - the 
Outer Lindblad Resonance - between the bar rotating at frequency Qb and the frequency of stars in the 
model Galactic disk. In particular, the training data that we use herein, was obtained by Chakrabarty 
(2007), by choosing the solar radial location from the interval [1.7, 2.3]. This explains the motivation for 
the aforementioned bounds on Rq, used in our work. 

To connect the used model to real astronomical units, we need to scale 0(i?) as in the simulation 
units, to the stellar rotational frequency, as measured in real astronomical units. This calls for the scaling 
of the constant vq to astronomical units, and of the galactocentric radial distance R, measured on this 
scale in which Rcr = 1, to astronomical units. The Sun is a star for which the galactocentric radial 
location can be learnt in units of Rcr using our inverse Bayesian learning strategy presented above. 
Now, the constant vq is known in real units as 220 kms^^ (Binney and Merrifield, 1998). Independent 
astronomical estimates of the galactocentric solar radial location have led to the standard value of 8 
kiloparsecs (Binney and Merrifield, 1998), where kiloparsec is an astronomical unit of distance. Then 
all distances encountered in the used model is connected to real astronomical unit of kiloparsec by 
scaling by the factor of 8/Rq kiloparsec, where Rq is the solar radial location that we learn here. 
The co-rotation radius is then scaled from the simulation units to real units as 1 x 8/Rq kiloparsecs. 
Then from the definition of co-rotation, the bar rotational frequency (If,, scaled to real units, stands at 
22O/(8/i20) kms~^kpc^^. Learning the rotational frequency of the bar is the ulterior benefit of learning 
the solar radial location as in our approach. 

The azimuthal separation between the long axis of the bar and the line that joins the Sun to the 
Galactic centre, is also suggested in past astronomical modelling work to be an acute angle (Chakrabarty, 
2007; Englmaier and Gerhard, 1999; Fux, 2001). Indeed, the training data used here was generated in 
simulations performed by Chakrabarty (2007), in which <1>0 is chosen from the interval [0, 90°]. This 
motivates the consideration of the interval of [0, 90°] for the azimuthal location of the Sun. 

Given the bounds on Rq and $q presented above, in our TMCMC algorithm, we reject those moves 
that suggest Rq and ^q values that fall outside these presented intervals. 

The 4 astrophysical models of the Galaxy that we discuss here are marked by the same choice of 
the value of Q^ and the background Galactic model parameters, while they are distinguished by the 
varying choices of the ratio $7^ : Qh, where the Galactic spiral pattern rotates with frequency 0^. In 
fact, the astrophysical model barQ is the only one that does not include the influence of the spiral pat- 
tern while the other three astrophysical models include the influence of both the bar and the spiral. 
In fact, the astrophysical models sp3bar3A8, sp3bar3 and sp3bar3J25 suggest Q,s '■ ^b values of 
18^1)/ 55, 22Qi)/55, 250;,/55. The physical effect of this choice is to induce varying levels of chaoticity 
in the 4 astrophysical models. Thus, Chakr^abarty and Sideris (2008) confirmed that of the 4 models, 
bar6 manifests very low chaoticity while sp3bar3 manifests maximal chaos at high energies, though 
both sp3bar3A8, sp3bar3J25 are comparably chaotic. 

Ancillary real data needs to be brought in to judge the relative fit amongst the astrophysical base 
models. In fact, Chakrabarty (2007) brought in extra information to perform model selection. Such in- 
formation was about the observed variance of the radial and azimuthal components of stellar velocities 
and this was used to rule out the model bar^ as physically viable, though the other three models were 
all acceptable from the point of view of all such ancillary observations that are available. This led to the 
inference that Og G [18^2^/55, 2517^/55]. It is to be noted that the 4 data sets simulated at the 4 dis- 
tinct astrophysical models, cannot be meaningfully combined since the data sets are based on mutually 
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2.2 2.3 
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Fig 1. Posteriors of Rq in units of RcR and $0 (in degrees) for the model bar6. 




Fig 2. Posteriors of Rq in units of RcR and $0 (in degrees) for the model spSbarS. 



contradicting base models of the system. Hence we enumerate the results from each of the 4 models 
individually. 

The marginal posterior densities of {Rq,^q) corresponding to the 4 astrophysical models of the 
Milky Way, are shown in Figures 1, 2, 3 and 4. Table 1 presents the posterior mode, the 95% highest 
posterior density (HPD) credible region of Rq and <I>q respectively, associated with the four models. 
Here Rq is expressed in the model units of length, i.e. in units of Rcr- ^q is expressed in degrees. 
The HPDs are computed using the methodology discussed in Carlin and Louis (1996). Disjoint HPD 
regions, characterise the highly multi-modal posterior distributions of the unknown location. Using the 
95% HPDs of the learnt Rq expressed in model units, and using the independently known astronomical 
measurement of the solar radial location as 8 kiloparsecs, the bar rotational frequency Q^ is computed 
(see explanation above) in Table 1. 

Summaries of the posteriors (mean, variance and 95% credible interval) of the smoothness parameters 
61, 62 and "E are presented in Tables 2, 3. Notable in all these tables are the small posterior variances 
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Fig 3. Posteriors of Rq in units of Rcr and <I>q (in degrees) for the model sp36ar3_18. 
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Fig 4. Posteriors of Rq in units of Rcr and $0 (in degrees) for the model sp3fear3_25. 



Table 1 
Summary of the posterior distributions of the radial component Rq and azimuthal component $0 of the unknown observer 
location vector for the 4 base astrophysical models and the unknown bar rotational frequency Q,b computed using the 95% 

HPDs on the learnt radial location Rq in these models. 



Model 


Rq (in units of Rcr) 


Qb (in kms ^/kiloparsec) 


$0 




Mode 


95% HPD 


95% HPD 


Mode 


95% HPD 


bar6 

spSbarS 

sp3bar 3_18 

sp36ar3.25 


2.20 
1.73 
1.76 
1.95 


[2.04, 2.30] 

[1.70,2.26] U [2.27,2.28] 

[1.70,2.29] 

[1.70,2.15] 


[56.1,63.25] 

[46.75, 62.15] U [62.45, 62.7] 

[46.75, 62.98] 

[46.75,59.12] 


23.50 
18.8 
32.5 
37.6 


[21.20,25.80] 

[9.6,61.5] 
[17.60, 79.90] 
[28.80, 40.40] 
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Table 2 
Summary of the posterior distributions of the smoothness parameters bi, b2for the 4 models. 



Model 


bi 


b2 




Mean Var 95% CI 


Mean Var 95% CI 


bar6 

spSbarS 

sp36ar3_18 

sp3fear3_25 


0.9598155 3.15x10"'' [0.959703,0.959879] 
0.8739616 6.72 x 10"^ [0.872347, 0.875052] 
0.9410686 1.46 x lO"'' [0.938852, 0.955264] 
0.7597931 5.64x10""' [0.759743,0.759833] 


1.005078 2.85 x 10"" [1.004985, 1.005142] 
1.003729 8.98 x 10"^ [1.002500, 1.005500] 
0.999010 4.08 x 10"® [0.997219, 1.004945] 
0.992174 2.89 x 10"^ [0.992067, 0.992246] 



Table 3 
Summary of the posterior distribution of the diagonal and one non-diagonal element ofS, from the 4 base astrophysical 

models. 



Model 


o-ii 


0"22 


CI"12 




95% CI 


95% CI 


95% CI 


bar6 

spSbarS 

sp3bar 3_18 

sp3bar 3_25 


[5.40 X 10"\4.0 X 10"*] 
[3.66 X 10"^1.03 X 10"^] 
[1.45 X 10"^1.68 X 10"^] 
[1.21 X 10"*, 5.69 X 10"*] 


[6.20 X 10"", 4.76 X 10"*] 
[6.53 X 10"^1.83 X 10"^] 
[1.29 X 10"^1.50 X 10"^] 
[1.13 X 10"*, 5.21 X 10"*] 


[0, 1.30 X 10""] 
[-6.40 X 10"", 2.68 X 10"*] 
[-1.19 X 10"*, 2.16 X 10"^] 
[-1.00 X 10"®, 1.50 X lO"'^] 



of the quantities in question; this is indicative of the fact that the data sets we used, in spite of the 
relatively smaller size compared to the astronomically large data sets used in the previous approaches in 
the literature, are very much informative, given our vector-variate QV-ba.sed Bayesian approach. 

As discussed in Section 5.6, thanks to our Gaussian Process approach, the posterior of XI should be 
close to the null matrix a posteriori if the choice of the design set and the number of design points are 
adequate. Quite encouragingly. Table 3, shows that indeed S is close to the null matrix a posteriori, for 
all the four models, signifying that the unknown velocity function has been learned well in all the cases. 

5.5. Comparison with results in astrophysical literature. The posterior summaries that we have pre- 
sented in the last section, for each of the 4 astrophysical models of the Milky Way, compare favourably 
with results obtained by Chakrabarty (2007), using the calibration method. However, the all important 
difference in our implementation is the vastly smaller data set that we needed to invoke, in order to 
achieve the learning of the two-dimensional vector S -in fact while in the calibration approach, the re- 
quired sample size is of the order of 3,500, in our work, this number is 50. Thus, data sufficiency issues, 
when a concern, are well tackled by our method. 

Upon the analyses of the viable astrophysical models of the Galaxy, Chakrabarty (2007) reported 
the result that Rq G [1.9375, 2. 21]i?c'iJ while $© G [0°,30°], where these ranges correspond to the 
presented uncertainties on the estimates, which were however, rather unsatisfactorally achieved (see 
Section 4). The values of the components of S, learnt in our work, overlap well with these results. As 
mentioned in Section 1, learning Rq and $0 have crucial astronomical significance for the azimuthal 
separation of the Sun from the long axis of the bar and for the rotational frequency ilf, of the bar in 
the Milky Way. That the models sp36ar3_18, sp3bar3 and sp3bar3J25 are distinguished by distinct 
values of the ratios of the rotational frequencies of the spiral pattern to the bar in the Galaxy, the derived 
estimate for Qi, (Table 1), in turn suggests possible credible regions on the value of the frequency Qg of 
the Milky Way spiral. 

However, the concurrence of our results with the results reported in astrophysical literature goes be- 
yond just the summaries of the posteriors of the solar position vector; remarkable con^elation can be no- 
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ticed between the measure of chaos in these 4 astrophysical models - as estimated by (Chakrabarty and Sideris, 
2008) - and the multi-modality of the posterior distribution of S that we advance. Chakrabarty and Sideris 
(2008) suggest minimum chaos in the bar6 model compared to the other three, while we notice the pos- 
teriors of both Rq and ^q in this model to be the unimodal - in fact the posteriors of Rq and $0 are 
unimodal only for this model, out of the 4 astrophysical models that we use to illustrate the efficacy 
of our method. Perhaps more importantly, the sp3bar3 model is noticed to manifest maximum (even 
global) chaoticity, on theoretical grounds by Chakrabarty (2007), backed by the chaos quantification at 
higher energies (Chakrabarty and Sideris, 2008). Likewise, in our work, the posterior distributions for 
Rq and ^0 are most multi-modal in this model, compared to the other three. The models sp3bar3A8 
and sp3bar3J25 are considered to be of intennediate chaoticity and we find these to correspond to pos- 
terior distributions (of sC"'^'")) that are multi-modal, though less so, than that for the model sp3bar3. 

The equivalent in the work by Chakrabarty (2007), of our marginal posterior distributions of {Rq, ^q)'^ 
is the multi-modality of the distribution of the p- value of the chosen test statistic in the space of S, as 
used by Chakrabarty (2007) to test for the null that the observed velocity data is drawn from the den- 
sity in the velocity space, estimated from the i-th simulated velocity data, for i = 1, . . . , n. When this 
distribution of the p-value was sought in the space of S, the distribution was indeed found to be most 
multi-modal for the sp3bar3 case while the modes were found to be less scattered than for the other 
three models. In line with this result, in our work we notice that the posterior probability of S is most 
multi-modal given the training data obtained with astrophysical model sp3bar3. The exact physical rea- 
son for the correlation between chaos in the astrophysical models and the multi-modality of the posterior 
distribution of S can be qualitatively understood if we begin with the premise that increased chaos gives 
rise to a more clustered distribution of the 2-dimensional stellar velocity vectors. This in turn implies 
that to an observer seated at any location on the Milky Way disk, the distribution of the velocities of n 
stars appears to be more multi-modal in nature than when there was less chaos (in the Galactic phase 
space). In other words, the distribution of the jTc-dimensional Vj vectors is more multi-modal when the 
training data set is constmcted using the sp3bar3 model than the bar_6 one. Consequently, the multi- 
modality in the achieved posterior probability density of s^"'^"') is higher in the former case than in the 
latter. 

Another point that merits mentions is that the estimates of Rq and $0 presented by Chakrabarty 
(2007) exclude the model sp3bar3 which could not be used to yield estimates given the highly scattered 
nature of the conesponding p- value distribution. Likewise, in our work, the same model manifests max- 
imal multi-modality amongst the others, but importantly, our approach allows for the representation of 
the full posterior density using which, the computation of the 95% HPDs is performed. 

5.6. Issue related to data sufficiency. It was noted in Section 5.4 above that in spite of our relatively 
smaller data sets, our learning of the unknown location vector is in accordance with the relevant results 
obtained by other approaches existing in the astrophysics literature. In this connection it is important 
to clarify how inference of the unknown quantities might be affected by the dimensionality of ^(•) 
which is in turn set by the number of stars chosen at any si, i = 1, . . . , n. Note that, iixespective of 
dimensionality, the Bayesian-ly predicted function must pass through all the points of the training data 
set. It is the I] matrix that bears information about covariance for a given star, (at a given value of the 
model parameter vector S), in our model. That the elements of the learnt S matrix are learnt to be small, 
and our learnt s("'^'^) are astrophysically interpretable, coiToborate the stability of our inversion given 
the n repeated samples. Such instabilities in inversion (Banerjee et al., 2008) can become effective for 
the following reason. If a particular functional value is supported by some argument, inclusion of more 
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(different) functions can change this value such that now it supports the old functions as well as the new 
functions. It is because we seek verification under such a situation that it becomes useful to learn the 
covariance matrix XI in the data (while the column covariance matrix C can be analytically integrated 
over). In fact, mai^ginalisation over C gives rise to a pai^ameter space that is more infoiTned by n than j, 
i.e. more informed by the number of known values of S at which data is simulated than the number of 
stars used at a given value of S. 

Also, each element of ^(•) has the same smoothness parameters for two given values of the Milky Way 
model pai-ameter vector (see Equation 4.9). Thus, increasing the number of stars while maintaining the 
number of design points a constant, will fail to improve inference about the smoothness parameters. We 
realise this by recalling that since the smoothness parameters determine the smoothness of the function 
^(•) with respect to the sample path - design points being the discretised version - information about the 
smoothness parameters is enhanced by increasing the number of design points rather than the number 
of stars. 

The arguments presented above show that the aim of improving the quality of inference on S and 
Q, is better served, if we increase the number n of design points in the used training data set, than the 
number of stars at each choice of the value Sj of the model parameter, i = 1, . . . , n. Understanding 
these data-sufficiency issues turns out to be very rewarding from the computational point of view - 
increasing the number of stars implies increasing the dimension of the matrix C, which would have 
made our computation prohibitively slow, but keeping the dimension of C at some reasonable level, and 
increasing n is computationally much less expensive. 

In our case, results for the unknown location are obtained by our approach with 50 x 2 velocity 
functions, which is a reasonable dimensionality, along with 216 well-chosen design points. Our results, 
which are astrophysically interpretable, demonstrate that these choices are adequate. 

5.7. Benefits of the method. An important benefit of our work is that it allows for the learning of 
a high-dimensional model parameter vector. Till now, a rigorous statistical learning of the Milky Way 
model parameter vector, from the inversion of available - real and simulated - stellai^ velocity data, has 
not taken place. Additionally, all previous attempts have acknowledged the need for including other 
Milky Way model parameters in their respective modelling strategies, but the lack of rigour in the same 
did not allow for the learning exercise to be performed in the context of a high-dimensional model 
parameter vector. For example, it is of substantial astronomical value if one could leam strengths of 
the perturbations due to the Galactic features as well other structural parameters and parameters that 
determine their dynamical behaviour (such as the ratio Q.s '■ ^b)- In our method, extension to such a 
high-dimensional S is easily possible, as long as astrophysical simulations carried out at varying values 
of these parameters are available. In fact, it is the matrix Ax>^^g (•) discussed above, that will be affected 
by the increase in the value of d; the computational complexity involved in the computation of this 
matrix scales linearly with d. 

The implication of our demonstrated success with a smaller data set is profound. Firstly, our method 
is advanced as a template for the analysis of the stellar phase space data that is available for the Milky 
Way, with the aim of learning a high-dimensional Galactic parameter vector; by extending the scope of 
the dynamical simulations of the Galaxy, performed on different astrophysical models of the Milky Way, 
the Milky Way models will be better constrained. At the same time, the planned flight of the GAIA - a 
mission of the European Space Agency - later in 2012, is set to provide large sets of phase space data all 
over the Milky Way. Our method, in conjunction with astrophysical models, can allow for fast learning 
of local and global model parameters of the Galaxy. In fact, the greater strength of our method lies in 
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its proposed applicability to galaxies other than our own, for which, only small data sets (< 100) can be 
expected. The learning of global model parameters of such systems is possible with a method such as 
ours, once data generated from simulations of such systems are available. 

6. Discussions. We have presented a new method to allow for inverse Bayesian learning of model 
parameter vectors, given real and simulated data, using a vector-variate Gaussian Process that was mo- 
tivated as an alternative for a matrix- variate QV of corresponding dimensions. The biggest advantage of 
our method is that it can work even when data sufficiency might be a concern and when making infer- 
ence about a high-dimensional model parameter vector is infeasible in other methods. The learning is 
rendered possible, owing to the matrix- variate Gaussian Process that we introduce here and successfully 
demonstrate an application of in an astronomical context. 

We saw in Section 5.7 that computational complexity scales only linearly with the dimensionality of 
the unknown model parameter S. Thus, porting a training data comprised of n independent values of 
S, Si, i = 1, . . . , n, where Sj is a d-dimensional vector, d > 2, is not going to render the computa- 
tional times infeasible. This allows for the learning of high-dimensional model parameter vectors in our 
method. 

In contrast to the situation with increasing the dimensionality of the unknown model parameter, in- 
creasing the dimensionality of the measurable will but imply substantial increase in the run time, since 
the relevant computational complexity then scales non-linearly, as about 0{k^), (in addition to the cost 
of k square roots), where k is the dimensionality of the measurable variable P. This is because of the 
dimensionality of the aforementioned S matrix is A; x A;, and the inverse of this enters the computation 
of the posterior. Thus, for example, increasing the measurable from a 2-dimensional to 4-dimensional 
vector increases the run time 8-fold, which is a large jump in the required run time. However, for most 
applications, we envisage the expansion of the dimensionality of the unknown model parameter, i.e. d, 
rather than that of the measurable, i.e. k. Thus, the method is expected to yield results within acceptable 
time frames, for most practical applications. 

The other major benefit of our work is that it allows for organic learning of the smoothness parameters, 
rather than results being subject to ad hoc choices of the same. 

Inverse learning similar to our application, is essentially of wide applicability - for example, one 
might need to learn the required dosage of a drug in order to achieve a cure efficacy that is judged 
optimal in curing a disease in a sample of j patients. This can be done using available data from a set 
of trials (the training data). Thus, in this scenario, n number of trials might be performed, each with 
j patients, such that changes in k characteristics of each patient are monitored where these measured 
"characteristics" reflect the efficacy of the administered dosage of the medication in affecting k different 
symptoms of the ailment. Thus, the efficacy parameter E is a fc-dimensional vector and we have data 
on the value of E of j patients, from n different trials that are distinguished by the n values of the 
dosage vector D. While dosage of the medication could be perceived as a scalar, it is also possible 
that it is rendered vector-valued, as can be the case, if there is reason to vary the dosage of the drug 
administered at different times of the day. In general, we consider D to be a d-dimensional vector. The 
test data is the data on the optimal improvement (cure efficacy) of k symptoms of these j patients when 
the administered dosage is the unknown Z)(*^'^*). Using the j x A;-dimensional efficacy matrices, attained 
for each of the n trials for the n values of the dosage vector, {Di}f^-^, Z)(*^'^*) can be learnt. Here, the 
matrix-valued information on the continuous parameter E (cure efficacy) is perceived as an unknown 
function of the dosage parameter vector which too is continuous in nature. We can realise this unknown 
function as a realisation of a Gaussian Process, after the method that we present above. 
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Alternatively, the price of a commodity might need to be fixed, given a parametrisation of the current 
demand for the item, guided by training data comprising such a demand parameter and the price, at time 
points in the past. 

In another astrophysical context, Chaki^abarty and Jackson (2009) learnt the distribution of the density 
of the total (i.e. luminous+dark) gravitational matter p : M.^ — > M, p{x, y, z) > 0, in distant galaxies 
belonging to an ubiquitous astronomical class of real galactic systems that is parametrised by the shape 
parameter n > and a length scale parameter Xe > 0. The measurements include the observed galactic 
image and central velocity dispersion measurement in the galaxy, while training data comprising syn- 
thetic image data and central velocity dispersion, generated from simulations of such systems earned out 
at chosen values of n, X^. and other relevant model variables. In such applications, as the measurables 
are rendered vectors and the data a matrix, inverse learning is possible using the method that we present 
here. 

As more Galactic simulations spanning a greater range of model parameters become available, the 
rigorous learning of such Milky Way parameters using our method will become possible, given the 
available stellar velocity data. This will enhance the quality of our knowledge about our own galaxy. That 
our method allows for such learning even for under-abundant systems, is encouraging for application of 
a similar analysis to galaxies other than our own, in which system parameters may be learnt using the 
much smaller available velocity data sets, compared to the situation in our galaxy. 

Supplementary material 

Section S-1 discusses the details of the dynamical simulations that lead to the training data set used in 
our supervised learning of the Milky Way model parameters. The posteriors of the matrices B and C 
that chai^acterise the Gaussian Process under consideration are presented in S-2. S-3 contains details of a 
proof that the matrix (n — m)kCGLS is positive definite. In Section S-4, a discussion about the batchwise 
updating within the TMCMC procedure is included. 
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SUPPLEMENTARY SECTION FOR "BAYESIAN NONPARAMETRIC ESTIMATION OF 
MILKY WAY MODEL PARAMETERS USING A NEW MATRIX- VARIATE GAUSSIAN 

PROCESS BASED METHOD" 

Throughout, we refer to our main manuscript as CBB. 

7, Details of dynamical simulations of astrophysical models. In Chakrabarty (2007), the mod- 
elling involves the following. A sample of phase space coordinates {wq}, is drawn from a chosen (to 
mimic real disc galaxies') phase space density g{w) at t = 0, and is evolved in a (chosen) paramet- 
ric Galactic gravitational potential ^ : R'^ x T — > M, where the time variable T G T C M>o. 
^(X, t) is chosen to emulate a realistic background Galactic potential ^o(x) € M, perturbed by 
the chosen parametric forms (rigidly rotating, quadrupolar) of the gravitational potential of the bar 
(^^(x) G M) and the (logarithmic spiral) gravitational potential of the spiral pattern (^^(x) G M) 
in the Milky Way. The perturbation strengths of these features (ef, : T — > M and e^ : T — > M) 
are slowly increased, to saturation values over time T^ G M, T^ > which is chosen such that 
Ts <^ Tsim, where Tgim G M is the total simulation time. Thus, in the q^^ astrophysical model, 
^(«)(x,t) = *[,''^(x) + eb{t)-^f\^,t) + es{m^^\^,t), for each g = 1, . . . ,4. At T = T,™,, or- 
bits are sampled in phase and recorded stroboscopically in the rotating frame of the bar at times when 
(^fe - O,)t=0. 

In order to sort the orbits by the value of the solar location in a galactocentric coordinate frame, 
chosen ranges of i? G [1.7,2.3] in simulation units and <1> G [0,90] in degrees, are used to define a 
regular, 2-D polar grid with Nr x N^ cells of radial and azimuthal widths of 5r and 5^, 5r, 6^ G M. 
The i*^' cell in this grid represents the i^^ value of the model parameter vector, S('\ with n = Nf> x Nff,, 
i = 1, . . . , Nf( X A'^^. In order to find stars that end up in the simulations with spatial coordinates within 
the interval around s^*) that defines the i*^-cell, we identify stars with radius G [s^ ,8^ + 6i) and 

azimuth G [s2 , S2 + <^2); these stars comprise the simulated velocity information Vj in the i*^-cell. 
This is performed for each i, i = 1, . . . , n. The same 2-D polar grid is used for each of the Milky Way 
astrophysical models. 

8. Posteriors of B and C given training data. Here, we discuss the achievement of the posteriors 
of the process matrices B and C, given the training data. To write the posteriors, we use uniform prior 
on B and a simple non-informative prior on C, namely, vr(C) oc| C |~(J+i)/2. 

Recalling that the distribution of the training data Vg is matrix normal, we write the posterior of B as 

(8.1) [B I i:,C,Q,Vs] r^ MMm,MBGLS,{HlA^^H^',n), 

where, Bqls '■= {hJ)A]j^Ho)'^{H^A~^^Vs); note that H^Bgls can be interpreted as the gener- 
alised least square (GLS) solution to the equation that given the chosen design vectors, the mean matrix 
HdB is equal to Vg. 

To write the posterior of C, for a chosen design set, we first define 

(8.2) Md := A^i - A^^Hd[hIA^^Hd]-^hIA^\ 
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and iV'^ M dVs)^^^^^^^ := [M*^^ i,t = 1, . . . , A:], where M*n is a matrix with j rows and j columns. 
Given Xl, we define m = d + 1 and V'j^^ as the {i, t)-th element of Xl^^, so that (n — m)kCGLs '■= 
Ylii=i St=i V'it ^*ti- Then it can be shown that for this given choice of the design vectors, the posterior 
distribution [C|5], Q^Vg] is inverse Wishart with parameters (n — m)kCGLS and {n — m)k (degrees 
of freedom). 

9. Positive deflniteness of (n — m)kCGLS- Now marginalisation of the posterior 

|-g(test) ^ ^^ Q^ -g^ C| v^*'^''*) , Vs] are sensitive to the posterior [C|S, Q, Vg]- Thus, we need to ensure that 
(n — rajkCcLS is positive definite in order for the posterior of C to be well-defined. We now present 
the proof that (n — m)kCGLS is positive definite. 

Proof. Since V^MV^ = [M*f.;i,t = l,...,k] is positive definite, B'^MV^ = PP', where 

P is a lower triangular matrix with strictly positive diagonal entries. Writing P = (P^:Pj^: • • • P'l)'^, 
where each Pt is of order j x j, it then follows that M*^ = P'^Pt- Hence, for any non-null j-component 
vector s, 

(9.1) J2E-u'^^M:,s = Y,Y.-r^Hp^snPts). 

i=l t=l i=l t=l 

Here an is the it-th element of the matrix (T£){-) := [a{-, si), . . . ,a{-, Sn)]^ The right hand side of the 
last equation is y'^{'S~^ Ij)y, where y := {yj, . . . , y^)^, with yf = Pts, and Ij is the j-th order 
identity matrix. Since Kronecker product of positive definite matrices yields positive definite matrices, 
it follows that the last expression is positive, implying that (n — m)kCGLS is positive definite. D 

10. Batchwise TMCMC based on additive transformations. Using TMCMC we can update the 
parameters in batches as follows: 

(i) Initialise the unknown quantities by fixing arbitrarily initial values is^^^^*'^\ Q^^' , S'^'^M. In our 

case, s^*'^***''^) = {si ' \ . . . , s|^ ^'^ ' '), Q^^' is chai^acterised by the initial values of the d smooth- 
ness parameters, which we denote by {b[ , . . . ,b^ ) and I]^*^' denotes the initial choice of the 
k X k matrix S. 
(ii) Assume that at iteration t, the state of the unknown parameters is (s(*est,t)^ Q^*s 5]^*^). 
(iii) (a) Propose e ~ 5'(-)^{e>o}' where g{-)is some arbitrary distribution, and / denotes the indicator 
function. In our applications, we shall choose g{-) = N{0, 1), so that, e > is drawn from 
a truncated normal distribution. Update the d components of location s^**^**) by setting, with 
probabilities vtj and (1 — ttj), s^^* ' = s-'^'^+Cje (forward transformation) and 

s^ "^^ ' ^ = s^- ^* ' — CjC (backward transformation), respectively, where, for j = 1, . . . , d, 
7Tj are appropriately chosen probabilities and Cj are appropriately chosen scaling factors. 
Assume that for ji G S, s^ '^^ ' ' gets the positive transformation, while for J2 G S'^, s^ '^* ' ' 
gets the backward transformation. Here S U S'^ = {I, . . . ,d}. 

(b) Accept the new proposal s(*^**'*+^) with acceptance probability 

(10.1) a^ = mm<l,— ^ — f- ^r x rA 
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where r^ denotes the ratio of An „„ ^Iff^ An -H'd^,,^ 2x51 2 (n+ 

A I J^aug I I JJaug Uaug ^aug I II IV' 

(71 + 1 — 77l)fc ^, //_l 1^ 

1 — m)kCGLS,aug\ 2 , evaluated at the new value (sl*est.*+iJ) and the current value 

^g(test,t)^ of s(*^**) respectively, with b and Xl remaining fixed at the cuixent values b^^' and 
I]^*\ respectively. In our applications, we shall choose ttj = 1/2 Vj, which simplifies the 
acceptance ratio by eliminating the first factor altogether, 
(c) If the new proposal s(*e^*'*+i) is not accepted, we set s(*e**'*+i) = g{test,t) _ 
(iv) Using the same TMCMC strategy, update b. Note that the probability of acceptance is zero if 

some components of the proposed value of b is negative. 
(v) Update S by first decomposing it into LL', where L is the appropriate lower-triangular matrix. 
Then update all the non-zero elements of Z in a single block using the same TMCMC strategy. The 
form of the acceptance ratio remains same. Again, reject the moves that yield negative diagonal 
elements. 
(vi) Cycle over steps (ii)-(v) for Ki + K2 iterations, assuming that the Markov chain converges after 
Ki iterations. Store the realisations 
{(«(*-*'*), Q(*),I]W) :t = Ki + l,...,K2] 

as samples obtained (approximately) from [s(*'^^*\Q,S | Vaug]- In particular; the realisations 
|g(test,t) . ^ _ j^j^ -)- 1^ . . . ,^2} are samples (approximately) from the marginal distribution 
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